Imaging a Black Hole using only Closure Quantities

In this tutorial, we will create a preliminary reconstruction of the 2017 M87 data on April 6 using closure-only imaging. This tutorial is a general introduction to closure-only imaging in Comrade. For an introduction to simultaneous image and instrument modeling, see Stokes I Simultaneous Image and Instrument Modeling

Introduction to Closure Imaging

The EHT is the highest-resolution telescope ever created. Its resolution is equivalent to roughly tracking a hockey puck on the moon when viewing it from the earth. However, the EHT is also a unique interferometer. For one, the data it produces is incredibly sparse. The array is formed from only eight geographic locations around the planet, each with its unique telescope. Additionally, the EHT observes at a much higher frequency than typical interferometers. As a result, it is often difficult to directly provide calibrated data since the source model can be complicated. This implies there can be large instrumental effects often called gains that can corrupt our signal. One way to deal with this is to fit quantities that are independent of gains. These are often called closure quantities. The types of closure quantities are briefly described in Introduction to the VLBI Imaging Problem.

In this tutorial, we will do closure-only modeling of M87 to produce preliminary images of M87.

To get started, we will load Comrade

using Comrade




using Pyehtim

  Activating project at `~/work/Comrade.jl/Comrade.jl/examples`

For reproducibility we use a stable random number genreator

using StableRNGs
rng = StableRNG(123)

StableRNGs.LehmerRNG(state=0x000000000000000000000000000000f7)

Load the Data

To download the data visit https://doi.org/10.25739/g85n-f134 To load the eht-imaging obsdata object we do:

obs = ehtim.obsdata.load_uvfits(joinpath(dirname(pathof(Comrade)), "..", "examples", "SR1_M87_2017_096_lo_hops_netcal_StokesI.uvfits"))

Python: <ehtim.obsdata.Obsdata object at 0x7fc68fb568f0>

Now we do some minor preprocessing:

• Scan average the data since the data have been preprocessed so that the gain phases are coherent.
• Add 1% systematic noise to deal with calibration issues that cause 1% non-closing errors.

obs = scan_average(obs).add_fractional_noise(0.015)

Python: <ehtim.obsdata.Obsdata object at 0x7fc68fb56e60>

Now, we extract our closure quantities from the EHT data set.

dlcamp, dcphase  = extract_table(obs, LogClosureAmplitudes(;snrcut=3), ClosurePhases(;snrcut=3))

(EHTObservation{Float64,Comrade.EHTLogClosureAmplitudeDatum{Float64}, ...}
  source: M87
  mjd: 57849
  frequency: 2.27070703125e11
  bandwidth: 1.856e9
  stations: [:AA, :AP, :AZ, :JC, :LM, :PV, :SM]
  nsamples: 128
, EHTObservation{Float64,Comrade.EHTClosurePhaseDatum{Float64}, ...}
  source: M87
  mjd: 57849
  frequency: 2.27070703125e11
  bandwidth: 1.856e9
  stations: [:AA, :AP, :AZ, :JC, :LM, :PV, :SM]
  nsamples: 152
)

Build the Model/Posterior

For our model, we will be using an image model that consists of a raster of point sources, convolved with some pulse or kernel to make a ContinuousImage object with it Comrade's. generic image model.

function sky(θ, metadata)
    (;fg, c, σimg) = θ
    (;K, meanpr, grid, cache) = metadata
    # Construct the image model we fix the flux to 0.6 Jy in this case
    cp = meanpr .+ σimg.*c.params
    rast = ((1-fg))*K(to_simplex(CenteredLR(), cp))
    img = IntensityMap(rast, grid)
    m = ContinuousImage(img, cache)
    # Add a large-scale gaussian to deal with the over-resolved mas flux
    g = modify(Gaussian(), Stretch(μas2rad(250.0), μas2rad(250.0)), Renormalize(fg))
    return m + g
end

sky (generic function with 1 method)

Now, let's set up our image model. The EHT's nominal resolution is 20-25 μas. Additionally, the EHT is not very sensitive to a larger field of views; typically, 60-80 μas is enough to describe the compact flux of M87. Given this, we only need to use a small number of pixels to describe our image.

npix = 32
fovxy = μas2rad(150.0)

7.27220521664304e-10

Now, we can feed in the array information to form the cache

grid = imagepixels(fovxy, fovxy, npix, npix)
buffer = IntensityMap(zeros(npix,npix), grid)
cache = create_cache(NFFTAlg(dlcamp), buffer, BSplinePulse{3}())

VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}(VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}(NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}(1, 4, 2.0, :kaiser_bessel, AbstractNFFTs.POLYNOMIAL, true, false, true, 0x00000000), [-4.405690154666661e9 787577.6145833326 … -5.999801315555549e9 -15551.297851562484; -4.523017159111106e9 -1.6838098888888871e6 … 3.059254300444441e9 118294.64453124987]), NFFTPlan with 274 sampling points for an input array of size(32, 32) and an output array of size(274,) with dims 1:2, ComplexF64[0.7014670924816703 - 0.519511453845776im, 0.9999999862106403 - 6.398624051637031e-5im, 0.6947985314160843 - 0.5233374146500747im, 0.9328845732113544 - 0.16336912604935844im, 0.9329421134238738 + 0.16331744240563614im, 0.6948519110447875 + 0.5233083088067162im, 0.7925669922789302 + 0.40478845348863063im, 0.9999999864389528 - 6.519719206334107e-5im, 0.9340555032742773 - 0.1520910238460368im, 0.6896791530194316 - 0.5261728587493705im  …  0.8376666586657548 + 0.1784956664688869im, 0.9196690775195503 + 0.09173305645579018im, 0.9135584686897926 + 0.29199716871629755im, 0.9135610039654034 + 0.2919905620278203im, 0.8899638273254515 + 0.28719490128071723im, 0.9568100109264257 + 0.09814170760184636im, 0.8899622135777404 + 0.287201492222944im, 0.9925668796456906 + 0.0027765725149323534im, 0.8376650833595861 - 0.17848899851834446im, 0.9999999999247222 + 7.335331210317892e-6im], BSplinePulse{3}(), [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0])

Now we need to specify our image prior. For this work we will use a Gaussian Markov Random field prior

using VLBIImagePriors, Distributions, DistributionsAD

Since we are using a Gaussian Markov random field prior we need to first specify our mean image. For this work we will use a symmetric Gaussian with a FWHM of 50 μas

fwhmfac = 2*sqrt(2*log(2))
mpr = modify(Gaussian(), Stretch(μas2rad(50.0)./fwhmfac))
imgpr = intensitymap(mpr, grid)

2-dimensional KeyedArray(NamedDimsArray(...)) with keys:
↓   X ∈ 32-element LinRange{Float64,...}
→   Y ∈ 32-element LinRange{Float64,...}
And data, 32×32 NamedDimsArray(::Matrix{Float64}, (:X, :Y)):
                 (-3.52247e-10)  …  (3.29522e-10)  (3.52247e-10)
 (-3.52247e-10)     6.37537e-8        1.32434e-7     6.37537e-8
 (-3.29522e-10)     1.32434e-7        2.751e-7       1.32434e-7
 (-3.06796e-10)     2.62014e-7        5.44273e-7     2.62014e-7
 (-2.84071e-10)     4.93724e-7        1.0256e-6      4.93724e-7
 (-2.61345e-10)     8.86092e-7   …    1.84065e-6     8.86092e-7
 (-2.38619e-10)     1.51463e-6        3.14629e-6     1.51463e-6
    ⋮                            ⋱    ⋮            
  (2.15894e-10)     2.46586e-6        5.12225e-6     2.46586e-6
  (2.38619e-10)     1.51463e-6        3.14629e-6     1.51463e-6
  (2.61345e-10)     8.86092e-7   …    1.84065e-6     8.86092e-7
  (2.84071e-10)     4.93724e-7        1.0256e-6      4.93724e-7
  (3.06796e-10)     2.62014e-7        5.44273e-7     2.62014e-7
  (3.29522e-10)     1.32434e-7        2.751e-7       1.32434e-7
  (3.52247e-10)     6.37537e-8        1.32434e-7     6.37537e-8

Now since we are actually modeling our image on the simplex we need to ensure that our mean image has unit flux

imgpr ./= flux(imgpr)

meanpr = to_real(CenteredLR(), Comrade.baseimage(imgpr))
metadata = (;meanpr,K=CenterImage(imgpr), grid, cache)

(meanpr = [-7.55422122563378 -6.823167558636962 … -6.823167558636962 -7.55422122563378; -6.823167558636962 -6.092113891640146 … -6.092113891640146 -6.823167558636962; … ; -6.823167558636962 -6.092113891640146 … -6.092113891640146 -6.823167558636962; -7.55422122563378 -6.823167558636962 … -6.823167558636962 -7.55422122563378], K = VLBIImagePriors.CenterImage{Matrix{Float64}, Tuple{Int64, Int64}}([0.9944957386363662 -0.005326704545454548 … 0.005326704545454572 0.005504261363636381; -0.005326704545454548 0.9948393969941349 … 0.005160603005865089 0.005326704545454565; … ; 0.005326704545454572 0.005160603005865089 … 0.9948393969941348 -0.005326704545454546; 0.005504261363636381 0.005326704545454565 … -0.005326704545454546 0.9944957386363638], (32, 32)), grid = GriddedKeys{(:X, :Y)}
	X: LinRange{Float64}(-3.5224744018114725e-10, 3.5224744018114725e-10, 32)
	Y: LinRange{Float64}(-3.5224744018114725e-10, 3.5224744018114725e-10, 32)
, cache = VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}(VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}(NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}(1, 4, 2.0, :kaiser_bessel, AbstractNFFTs.POLYNOMIAL, true, false, true, 0x00000000), [-4.405690154666661e9 787577.6145833326 … -5.999801315555549e9 -15551.297851562484; -4.523017159111106e9 -1.6838098888888871e6 … 3.059254300444441e9 118294.64453124987]), NFFTPlan with 274 sampling points for an input array of size(32, 32) and an output array of size(274,) with dims 1:2, ComplexF64[0.7014670924816703 - 0.519511453845776im, 0.9999999862106403 - 6.398624051637031e-5im, 0.6947985314160843 - 0.5233374146500747im, 0.9328845732113544 - 0.16336912604935844im, 0.9329421134238738 + 0.16331744240563614im, 0.6948519110447875 + 0.5233083088067162im, 0.7925669922789302 + 0.40478845348863063im, 0.9999999864389528 - 6.519719206334107e-5im, 0.9340555032742773 - 0.1520910238460368im, 0.6896791530194316 - 0.5261728587493705im  …  0.8376666586657548 + 0.1784956664688869im, 0.9196690775195503 + 0.09173305645579018im, 0.9135584686897926 + 0.29199716871629755im, 0.9135610039654034 + 0.2919905620278203im, 0.8899638273254515 + 0.28719490128071723im, 0.9568100109264257 + 0.09814170760184636im, 0.8899622135777404 + 0.287201492222944im, 0.9925668796456906 + 0.0027765725149323534im, 0.8376650833595861 - 0.17848899851834446im, 0.9999999999247222 + 7.335331210317892e-6im], BSplinePulse{3}(), [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]))

In addition we want a reasonable guess for what the resolution of our image should be. For radio astronomy this is given by roughly the longest baseline in the image. To put this into pixel space we then divide by the pixel size.

beam = beamsize(dlcamp)
rat = (beam/(step(grid.X)))

5.326336637737519

To make the Gaussian Markov random field efficient we first precompute a bunch of quantities that allow us to scale things linearly with the number of image pixels. This drastically improves the usual N^3 scaling you get from usual Gaussian Processes.

crcache = MarkovRandomFieldCache(meanpr)

VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}(sparse([1, 2, 32, 33, 993, 1, 2, 3, 34, 994  …  31, 991, 1022, 1023, 1024, 32, 992, 993, 1023, 1024], [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 13, 13, 13, 13, 13, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15, 16, 16, 16, 16, 16, 17, 17, 17, 17, 17, 18, 18, 18, 18, 18, 19, 19, 19, 19, 19, 20, 20, 20, 20, 20, 21, 21, 21, 21, 21, 22, 22, 22, 22, 22, 23, 23, 23, 23, 23, 24, 24, 24, 24, 24, 25, 25, 25, 25, 25, 26, 26, 26, 26, 26, 27, 27, 27, 27, 27, 28, 28, 28, 28, 28, 29, 29, 29, 29, 29, 30, 30, 30, 30, 30, 31, 31, 31, 31, 31, 32, 32, 32, 32, 32, 33, 33, 33, 33, 33, 34, 34, 34, 34, 34, 35, 35, 35, 35, 35, 36, 36, 36, 36, 36, 37, 37, 37, 37, 37, 38, 38, 38, 38, 38, 39, 39, 39, 39, 39, 40, 40, 40, 40, 40, 41, 41, 41, 41, 41, 42, 42, 42, 42, 42, 43, 43, 43, 43, 43, 44, 44, 44, 44, 44, 45, 45, 45, 45, 45, 46, 46, 46, 46, 46, 47, 47, 47, 47, 47, 48, 48, 48, 48, 48, 49, 49, 49, 49, 49, 50, 50, 50, 50, 50, 51, 51, 51, 51, 51, 52, 52, 52, 52, 52, 53, 53, 53, 53, 53, 54, 54, 54, 54, 54, 55, 55, 55, 55, 55, 56, 56, 56, 56, 56, 57, 57, 57, 57, 57, 58, 58, 58, 58, 58, 59, 59, 59, 59, 59, 60, 60, 60, 60, 60, 61, 61, 61, 61, 61, 62, 62, 62, 62, 62, 63, 63, 63, 63, 63, 64, 64, 64, 64, 64, 65, 65, 65, 65, 65, 66, 66, 66, 66, 66, 67, 67, 67, 67, 67, 68, 68, 68, 68, 68, 69, 69, 69, 69, 69, 70, 70, 70, 70, 70, 71, 71, 71, 71, 71, 72, 72, 72, 72, 72, 73, 73, 73, 73, 73, 74, 74, 74, 74, 74, 75, 75, 75, 75, 75, 76, 76, 76, 76, 76, 77, 77, 77, 77, 77, 78, 78, 78, 78, 78, 79, 79, 79, 79, 79, 80, 80, 80, 80, 80, 81, 81, 81, 81, 81, 82, 82, 82, 82, 82, 83, 83, 83, 83, 83, 84, 84, 84, 84, 84, 85, 85, 85, 85, 85, 86, 86, 86, 86, 86, 87, 87, 87, 87, 87, 88, 88, 88, 88, 88, 89, 89, 89, 89, 89, 90, 90, 90, 90, 90, 91, 91, 91, 91, 91, 92, 92, 92, 92, 92, 93, 93, 93, 93, 93, 94, 94, 94, 94, 94, 95, 95, 95, 95, 95, 96, 96, 96, 96, 96, 97, 97, 97, 97, 97, 98, 98, 98, 98, 98, 99, 99, 99, 99, 99, 100, 100, 100, 100, 100, 101, 101, 101, 101, 101, 102, 102, 102, 102, 102, 103, 103, 103, 103, 103, 104, 104, 104, 104, 104, 105, 105, 105, 105, 105, 106, 106, 106, 106, 106, 107, 107, 107, 107, 107, 108, 108, 108, 108, 108, 109, 109, 109, 109, 109, 110, 110, 110, 110, 110, 111, 111, 111, 111, 111, 112, 112, 112, 112, 112, 113, 113, 113, 113, 113, 114, 114, 114, 114, 114, 115, 115, 115, 115, 115, 116, 116, 116, 116, 116, 117, 117, 117, 117, 117, 118, 118, 118, 118, 118, 119, 119, 119, 119, 119, 120, 120, 120, 120, 120, 121, 121, 121, 121, 121, 122, 122, 122, 122, 122, 123, 123, 123, 123, 123, 124, 124, 124, 124, 124, 125, 125, 125, 125, 125, 126, 126, 126, 126, 126, 127, 127, 127, 127, 127, 128, 128, 128, 128, 128, 129, 129, 129, 129, 129, 130, 130, 130, 130, 130, 131, 131, 131, 131, 131, 132, 132, 132, 132, 132, 133, 133, 133, 133, 133, 134, 134, 134, 134, 134, 135, 135, 135, 135, 135, 136, 136, 136, 136, 136, 137, 137, 137, 137, 137, 138, 138, 138, 138, 138, 139, 139, 139, 139, 139, 140, 140, 140, 140, 140, 141, 141, 141, 141, 141, 142, 142, 142, 142, 142, 143, 143, 143, 143, 143, 144, 144, 144, 144, 144, 145, 145, 145, 145, 145, 146, 146, 146, 146, 146, 147, 147, 147, 147, 147, 148, 148, 148, 148, 148, 149, 149, 149, 149, 149, 150, 150, 150, 150, 150, 151, 151, 151, 151, 151, 152, 152, 152, 152, 152, 153, 153, 153, 153, 153, 154, 154, 154, 154, 154, 155, 155, 155, 155, 155, 156, 156, 156, 156, 156, 157, 157, 157, 157, 157, 158, 158, 158, 158, 158, 159, 159, 159, 159, 159, 160, 160, 160, 160, 160, 161, 161, 161, 161, 161, 162, 162, 162, 162, 162, 163, 163, 163, 163, 163, 164, 164, 164, 164, 164, 165, 165, 165, 165, 165, 166, 166, 166, 166, 166, 167, 167, 167, 167, 167, 168, 168, 168, 168, 168, 169, 169, 169, 169, 169, 170, 170, 170, 170, 170, 171, 171, 171, 171, 171, 172, 172, 172, 172, 172, 173, 173, 173, 173, 173, 174, 174, 174, 174, 174, 175, 175, 175, 175, 175, 176, 176, 176, 176, 176, 177, 177, 177, 177, 177, 178, 178, 178, 178, 178, 179, 179, 179, 179, 179, 180, 180, 180, 180, 180, 181, 181, 181, 181, 181, 182, 182, 182, 182, 182, 183, 183, 183, 183, 183, 184, 184, 184, 184, 184, 185, 185, 185, 185, 185, 186, 186, 186, 186, 186, 187, 187, 187, 187, 187, 188, 188, 188, 188, 188, 189, 189, 189, 189, 189, 190, 190, 190, 190, 190, 191, 191, 191, 191, 191, 192, 192, 192, 192, 192, 193, 193, 193, 193, 193, 194, 194, 194, 194, 194, 195, 195, 195, 195, 195, 196, 196, 196, 196, 196, 197, 197, 197, 197, 197, 198, 198, 198, 198, 198, 199, 199, 199, 199, 199, 200, 200, 200, 200, 200, 201, 201, 201, 201, 201, 202, 202, 202, 202, 202, 203, 203, 203, 203, 203, 204, 204, 204, 204, 204, 205, 205, 205, 205, 205, 206, 206, 206, 206, 206, 207, 207, 207, 207, 207, 208, 208, 208, 208, 208, 209, 209, 209, 209, 209, 210, 210, 210, 210, 210, 211, 211, 211, 211, 211, 212, 212, 212, 212, 212, 213, 213, 213, 213, 213, 214, 214, 214, 214, 214, 215, 215, 215, 215, 215, 216, 216, 216, 216, 216, 217, 217, 217, 217, 217, 218, 218, 218, 218, 218, 219, 219, 219, 219, 219, 220, 220, 220, 220, 220, 221, 221, 221, 221, 221, 222, 222, 222, 222, 222, 223, 223, 223, 223, 223, 224, 224, 224, 224, 224, 225, 225, 225, 225, 225, 226, 226, 226, 226, 226, 227, 227, 227, 227, 227, 228, 228, 228, 228, 228, 229, 229, 229, 229, 229, 230, 230, 230, 230, 230, 231, 231, 231, 231, 231, 232, 232, 232, 232, 232, 233, 233, 233, 233, 233, 234, 234, 234, 234, 234, 235, 235, 235, 235, 235, 236, 236, 236, 236, 236, 237, 237, 237, 237, 237, 238, 238, 238, 238, 238, 239, 239, 239, 239, 239, 240, 240, 240, 240, 240, 241, 241, 241, 241, 241, 242, 242, 242, 242, 242, 243, 243, 243, 243, 243, 244, 244, 244, 244, 244, 245, 245, 245, 245, 245, 246, 246, 246, 246, 246, 247, 247, 247, 247, 247, 248, 248, 248, 248, 248, 249, 249, 249, 249, 249, 250, 250, 250, 250, 250, 251, 251, 251, 251, 251, 252, 252, 252, 252, 252, 253, 253, 253, 253, 253, 254, 254, 254, 254, 254, 255, 255, 255, 255, 255, 256, 256, 256, 256, 256, 257, 257, 257, 257, 257, 258, 258, 258, 258, 258, 259, 259, 259, 259, 259, 260, 260, 260, 260, 260, 261, 261, 261, 261, 261, 262, 262, 262, 262, 262, 263, 263, 263, 263, 263, 264, 264, 264, 264, 264, 265, 265, 265, 265, 265, 266, 266, 266, 266, 266, 267, 267, 267, 267, 267, 268, 268, 268, 268, 268, 269, 269, 269, 269, 269, 270, 270, 270, 270, 270, 271, 271, 271, 271, 271, 272, 272, 272, 272, 272, 273, 273, 273, 273, 273, 274, 274, 274, 274, 274, 275, 275, 275, 275, 275, 276, 276, 276, 276, 276, 277, 277, 277, 277, 277, 278, 278, 278, 278, 278, 279, 279, 279, 279, 279, 280, 280, 280, 280, 280, 281, 281, 281, 281, 281, 282, 282, 282, 282, 282, 283, 283, 283, 283, 283, 284, 284, 284, 284, 284, 285, 285, 285, 285, 285, 286, 286, 286, 286, 286, 287, 287, 287, 287, 287, 288, 288, 288, 288, 288, 289, 289, 289, 289, 289, 290, 290, 290, 290, 290, 291, 291, 291, 291, 291, 292, 292, 292, 292, 292, 293, 293, 293, 293, 293, 294, 294, 294, 294, 294, 295, 295, 295, 295, 295, 296, 296, 296, 296, 296, 297, 297, 297, 297, 297, 298, 298, 298, 298, 298, 299, 299, 299, 299, 299, 300, 300, 300, 300, 300, 301, 301, 301, 301, 301, 302, 302, 302, 302, 302, 303, 303, 303, 303, 303, 304, 304, 304, 304, 304, 305, 305, 305, 305, 305, 306, 306, 306, 306, 306, 307, 307, 307, 307, 307, 308, 308, 308, 308, 308, 309, 309, 309, 309, 309, 310, 310, 310, 310, 310, 311, 311, 311, 311, 311, 312, 312, 312, 312, 312, 313, 313, 313, 313, 313, 314, 314, 314, 314, 314, 315, 315, 315, 315, 315, 316, 316, 316, 316, 316, 317, 317, 317, 317, 317, 318, 318, 318, 318, 318, 319, 319, 319, 319, 319, 320, 320, 320, 320, 320, 321, 321, 321, 321, 321, 322, 322, 322, 322, 322, 323, 323, 323, 323, 323, 324, 324, 324, 324, 324, 325, 325, 325, 325, 325, 326, 326, 326, 326, 326, 327, 327, 327, 327, 327, 328, 328, 328, 328, 328, 329, 329, 329, 329, 329, 330, 330, 330, 330, 330, 331, 331, 331, 331, 331, 332, 332, 332, 332, 332, 333, 333, 333, 333, 333, 334, 334, 334, 334, 334, 335, 335, 335, 335, 335, 336, 336, 336, 336, 336, 337, 337, 337, 337, 337, 338, 338, 338, 338, 338, 339, 339, 339, 339, 339, 340, 340, 340, 340, 340, 341, 341, 341, 341, 341, 342, 342, 342, 342, 342, 343, 343, 343, 343, 343, 344, 344, 344, 344, 344, 345, 345, 345, 345, 345, 346, 346, 346, 346, 346, 347, 347, 347, 347, 347, 348, 348, 348, 348, 348, 349, 349, 349, 349, 349, 350, 350, 350, 350, 350, 351, 351, 351, 351, 351, 352, 352, 352, 352, 352, 353, 353, 353, 353, 353, 354, 354, 354, 354, 354, 355, 355, 355, 355, 355, 356, 356, 356, 356, 356, 357, 357, 357, 357, 357, 358, 358, 358, 358, 358, 359, 359, 359, 359, 359, 360, 360, 360, 360, 360, 361, 361, 361, 361, 361, 362, 362, 362, 362, 362, 363, 363, 363, 363, 363, 364, 364, 364, 364, 364, 365, 365, 365, 365, 365, 366, 366, 366, 366, 366, 367, 367, 367, 367, 367, 368, 368, 368, 368, 368, 369, 369, 369, 369, 369, 370, 370, 370, 370, 370, 371, 371, 371, 371, 371, 372, 372, 372, 372, 372, 373, 373, 373, 373, 373, 374, 374, 374, 374, 374, 375, 375, 375, 375, 375, 376, 376, 376, 376, 376, 377, 377, 377, 377, 377, 378, 378, 378, 378, 378, 379, 379, 379, 379, 379, 380, 380, 380, 380, 380, 381, 381, 381, 381, 381, 382, 382, 382, 382, 382, 383, 383, 383, 383, 383, 384, 384, 384, 384, 384, 385, 385, 385, 385, 385, 386, 386, 386, 386, 386, 387, 387, 387, 387, 387, 388, 388, 388, 388, 388, 389, 389, 389, 389, 389, 390, 390, 390, 390, 390, 391, 391, 391, 391, 391, 392, 392, 392, 392, 392, 393, 393, 393, 393, 393, 394, 394, 394, 394, 394, 395, 395, 395, 395, 395, 396, 396, 396, 396, 396, 397, 397, 397, 397, 397, 398, 398, 398, 398, 398, 399, 399, 399, 399, 399, 400, 400, 400, 400, 400, 401, 401, 401, 401, 401, 402, 402, 402, 402, 402, 403, 403, 403, 403, 403, 404, 404, 404, 404, 404, 405, 405, 405, 405, 405, 406, 406, 406, 406, 406, 407, 407, 407, 407, 407, 408, 408, 408, 408, 408, 409, 409, 409, 409, 409, 410, 410, 410, 410, 410, 411, 411, 411, 411, 411, 412, 412, 412, 412, 412, 413, 413, 413, 413, 413, 414, 414, 414, 414, 414, 415, 415, 415, 415, 415, 416, 416, 416, 416, 416, 417, 417, 417, 417, 417, 418, 418, 418, 418, 418, 419, 419, 419, 419, 419, 420, 420, 420, 420, 420, 421, 421, 421, 421, 421, 422, 422, 422, 422, 422, 423, 423, 423, 423, 423, 424, 424, 424, 424, 424, 425, 425, 425, 425, 425, 426, 426, 426, 426, 426, 427, 427, 427, 427, 427, 428, 428, 428, 428, 428, 429, 429, 429, 429, 429, 430, 430, 430, 430, 430, 431, 431, 431, 431, 431, 432, 432, 432, 432, 432, 433, 433, 433, 433, 433, 434, 434, 434, 434, 434, 435, 435, 435, 435, 435, 436, 436, 436, 436, 436, 437, 437, 437, 437, 437, 438, 438, 438, 438, 438, 439, 439, 439, 439, 439, 440, 440, 440, 440, 440, 441, 441, 441, 441, 441, 442, 442, 442, 442, 442, 443, 443, 443, 443, 443, 444, 444, 444, 444, 444, 445, 445, 445, 445, 445, 446, 446, 446, 446, 446, 447, 447, 447, 447, 447, 448, 448, 448, 448, 448, 449, 449, 449, 449, 449, 450, 450, 450, 450, 450, 451, 451, 451, 451, 451, 452, 452, 452, 452, 452, 453, 453, 453, 453, 453, 454, 454, 454, 454, 454, 455, 455, 455, 455, 455, 456, 456, 456, 456, 456, 457, 457, 457, 457, 457, 458, 458, 458, 458, 458, 459, 459, 459, 459, 459, 460, 460, 460, 460, 460, 461, 461, 461, 461, 461, 462, 462, 462, 462, 462, 463, 463, 463, 463, 463, 464, 464, 464, 464, 464, 465, 465, 465, 465, 465, 466, 466, 466, 466, 466, 467, 467, 467, 467, 467, 468, 468, 468, 468, 468, 469, 469, 469, 469, 469, 470, 470, 470, 470, 470, 471, 471, 471, 471, 471, 472, 472, 472, 472, 472, 473, 473, 473, 473, 473, 474, 474, 474, 474, 474, 475, 475, 475, 475, 475, 476, 476, 476, 476, 476, 477, 477, 477, 477, 477, 478, 478, 478, 478, 478, 479, 479, 479, 479, 479, 480, 480, 480, 480, 480, 481, 481, 481, 481, 481, 482, 482, 482, 482, 482, 483, 483, 483, 483, 483, 484, 484, 484, 484, 484, 485, 485, 485, 485, 485, 486, 486, 486, 486, 486, 487, 487, 487, 487, 487, 488, 488, 488, 488, 488, 489, 489, 489, 489, 489, 490, 490, 490, 490, 490, 491, 491, 491, 491, 491, 492, 492, 492, 492, 492, 493, 493, 493, 493, 493, 494, 494, 494, 494, 494, 495, 495, 495, 495, 495, 496, 496, 496, 496, 496, 497, 497, 497, 497, 497, 498, 498, 498, 498, 498, 499, 499, 499, 499, 499, 500, 500, 500, 500, 500, 501, 501, 501, 501, 501, 502, 502, 502, 502, 502, 503, 503, 503, 503, 503, 504, 504, 504, 504, 504, 505, 505, 505, 505, 505, 506, 506, 506, 506, 506, 507, 507, 507, 507, 507, 508, 508, 508, 508, 508, 509, 509, 509, 509, 509, 510, 510, 510, 510, 510, 511, 511, 511, 511, 511, 512, 512, 512, 512, 512, 513, 513, 513, 513, 513, 514, 514, 514, 514, 514, 515, 515, 515, 515, 515, 516, 516, 516, 516, 516, 517, 517, 517, 517, 517, 518, 518, 518, 518, 518, 519, 519, 519, 519, 519, 520, 520, 520, 520, 520, 521, 521, 521, 521, 521, 522, 522, 522, 522, 522, 523, 523, 523, 523, 523, 524, 524, 524, 524, 524, 525, 525, 525, 525, 525, 526, 526, 526, 526, 526, 527, 527, 527, 527, 527, 528, 528, 528, 528, 528, 529, 529, 529, 529, 529, 530, 530, 530, 530, 530, 531, 531, 531, 531, 531, 532, 532, 532, 532, 532, 533, 533, 533, 533, 533, 534, 534, 534, 534, 534, 535, 535, 535, 535, 535, 536, 536, 536, 536, 536, 537, 537, 537, 537, 537, 538, 538, 538, 538, 538, 539, 539, 539, 539, 539, 540, 540, 540, 540, 540, 541, 541, 541, 541, 541, 542, 542, 542, 542, 542, 543, 543, 543, 543, 543, 544, 544, 544, 544, 544, 545, 545, 545, 545, 545, 546, 546, 546, 546, 546, 547, 547, 547, 547, 547, 548, 548, 548, 548, 548, 549, 549, 549, 549, 549, 550, 550, 550, 550, 550, 551, 551, 551, 551, 551, 552, 552, 552, 552, 552, 553, 553, 553, 553, 553, 554, 554, 554, 554, 554, 555, 555, 555, 555, 555, 556, 556, 556, 556, 556, 557, 557, 557, 557, 557, 558, 558, 558, 558, 558, 559, 559, 559, 559, 559, 560, 560, 560, 560, 560, 561, 561, 561, 561, 561, 562, 562, 562, 562, 562, 563, 563, 563, 563, 563, 564, 564, 564, 564, 564, 565, 565, 565, 565, 565, 566, 566, 566, 566, 566, 567, 567, 567, 567, 567, 568, 568, 568, 568, 568, 569, 569, 569, 569, 569, 570, 570, 570, 570, 570, 571, 571, 571, 571, 571, 572, 572, 572, 572, 572, 573, 573, 573, 573, 573, 574, 574, 574, 574, 574, 575, 575, 575, 575, 575, 576, 576, 576, 576, 576, 577, 577, 577, 577, 577, 578, 578, 578, 578, 578, 579, 579, 579, 579, 579, 580, 580, 580, 580, 580, 581, 581, 581, 581, 581, 582, 582, 582, 582, 582, 583, 583, 583, 583, 583, 584, 584, 584, 584, 584, 585, 585, 585, 585, 585, 586, 586, 586, 586, 586, 587, 587, 587, 587, 587, 588, 588, 588, 588, 588, 589, 589, 589, 589, 589, 590, 590, 590, 590, 590, 591, 591, 591, 591, 591, 592, 592, 592, 592, 592, 593, 593, 593, 593, 593, 594, 594, 594, 594, 594, 595, 595, 595, 595, 595, 596, 596, 596, 596, 596, 597, 597, 597, 597, 597, 598, 598, 598, 598, 598, 599, 599, 599, 599, 599, 600, 600, 600, 600, 600, 601, 601, 601, 601, 601, 602, 602, 602, 602, 602, 603, 603, 603, 603, 603, 604, 604, 604, 604, 604, 605, 605, 605, 605, 605, 606, 606, 606, 606, 606, 607, 607, 607, 607, 607, 608, 608, 608, 608, 608, 609, 609, 609, 609, 609, 610, 610, 610, 610, 610, 611, 611, 611, 611, 611, 612, 612, 612, 612, 612, 613, 613, 613, 613, 613, 614, 614, 614, 614, 614, 615, 615, 615, 615, 615, 616, 616, 616, 616, 616, 617, 617, 617, 617, 617, 618, 618, 618, 618, 618, 619, 619, 619, 619, 619, 620, 620, 620, 620, 620, 621, 621, 621, 621, 621, 622, 622, 622, 622, 622, 623, 623, 623, 623, 623, 624, 624, 624, 624, 624, 625, 625, 625, 625, 625, 626, 626, 626, 626, 626, 627, 627, 627, 627, 627, 628, 628, 628, 628, 628, 629, 629, 629, 629, 629, 630, 630, 630, 630, 630, 631, 631, 631, 631, 631, 632, 632, 632, 632, 632, 633, 633, 633, 633, 633, 634, 634, 634, 634, 634, 635, 635, 635, 635, 635, 636, 636, 636, 636, 636, 637, 637, 637, 637, 637, 638, 638, 638, 638, 638, 639, 639, 639, 639, 639, 640, 640, 640, 640, 640, 641, 641, 641, 641, 641, 642, 642, 642, 642, 642, 643, 643, 643, 643, 643, 644, 644, 644, 644, 644, 645, 645, 645, 645, 645, 646, 646, 646, 646, 646, 647, 647, 647, 647, 647, 648, 648, 648, 648, 648, 649, 649, 649, 649, 649, 650, 650, 650, 650, 650, 651, 651, 651, 651, 651, 652, 652, 652, 652, 652, 653, 653, 653, 653, 653, 654, 654, 654, 654, 654, 655, 655, 655, 655, 655, 656, 656, 656, 656, 656, 657, 657, 657, 657, 657, 658, 658, 658, 658, 658, 659, 659, 659, 659, 659, 660, 660, 660, 660, 660, 661, 661, 661, 661, 661, 662, 662, 662, 662, 662, 663, 663, 663, 663, 663, 664, 664, 664, 664, 664, 665, 665, 665, 665, 665, 666, 666, 666, 666, 666, 667, 667, 667, 667, 667, 668, 668, 668, 668, 668, 669, 669, 669, 669, 669, 670, 670, 670, 670, 670, 671, 671, 671, 671, 671, 672, 672, 672, 672, 672, 673, 673, 673, 673, 673, 674, 674, 674, 674, 674, 675, 675, 675, 675, 675, 676, 676, 676, 676, 676, 677, 677, 677, 677, 677, 678, 678, 678, 678, 678, 679, 679, 679, 679, 679, 680, 680, 680, 680, 680, 681, 681, 681, 681, 681, 682, 682, 682, 682, 682, 683, 683, 683, 683, 683, 684, 684, 684, 684, 684, 685, 685, 685, 685, 685, 686, 686, 686, 686, 686, 687, 687, 687, 687, 687, 688, 688, 688, 688, 688, 689, 689, 689, 689, 689, 690, 690, 690, 690, 690, 691, 691, 691, 691, 691, 692, 692, 692, 692, 692, 693, 693, 693, 693, 693, 694, 694, 694, 694, 694, 695, 695, 695, 695, 695, 696, 696, 696, 696, 696, 697, 697, 697, 697, 697, 698, 698, 698, 698, 698, 699, 699, 699, 699, 699, 700, 700, 700, 700, 700, 701, 701, 701, 701, 701, 702, 702, 702, 702, 702, 703, 703, 703, 703, 703, 704, 704, 704, 704, 704, 705, 705, 705, 705, 705, 706, 706, 706, 706, 706, 707, 707, 707, 707, 707, 708, 708, 708, 708, 708, 709, 709, 709, 709, 709, 710, 710, 710, 710, 710, 711, 711, 711, 711, 711, 712, 712, 712, 712, 712, 713, 713, 713, 713, 713, 714, 714, 714, 714, 714, 715, 715, 715, 715, 715, 716, 716, 716, 716, 716, 717, 717, 717, 717, 717, 718, 718, 718, 718, 718, 719, 719, 719, 719, 719, 720, 720, 720, 720, 720, 721, 721, 721, 721, 721, 722, 722, 722, 722, 722, 723, 723, 723, 723, 723, 724, 724, 724, 724, 724, 725, 725, 725, 725, 725, 726, 726, 726, 726, 726, 727, 727, 727, 727, 727, 728, 728, 728, 728, 728, 729, 729, 729, 729, 729, 730, 730, 730, 730, 730, 731, 731, 731, 731, 731, 732, 732, 732, 732, 732, 733, 733, 733, 733, 733, 734, 734, 734, 734, 734, 735, 735, 735, 735, 735, 736, 736, 736, 736, 736, 737, 737, 737, 737, 737, 738, 738, 738, 738, 738, 739, 739, 739, 739, 739, 740, 740, 740, 740, 740, 741, 741, 741, 741, 741, 742, 742, 742, 742, 742, 743, 743, 743, 743, 743, 744, 744, 744, 744, 744, 745, 745, 745, 745, 745, 746, 746, 746, 746, 746, 747, 747, 747, 747, 747, 748, 748, 748, 748, 748, 749, 749, 749, 749, 749, 750, 750, 750, 750, 750, 751, 751, 751, 751, 751, 752, 752, 752, 752, 752, 753, 753, 753, 753, 753, 754, 754, 754, 754, 754, 755, 755, 755, 755, 755, 756, 756, 756, 756, 756, 757, 757, 757, 757, 757, 758, 758, 758, 758, 758, 759, 759, 759, 759, 759, 760, 760, 760, 760, 760, 761, 761, 761, 761, 761, 762, 762, 762, 762, 762, 763, 763, 763, 763, 763, 764, 764, 764, 764, 764, 765, 765, 765, 765, 765, 766, 766, 766, 766, 766, 767, 767, 767, 767, 767, 768, 768, 768, 768, 768, 769, 769, 769, 769, 769, 770, 770, 770, 770, 770, 771, 771, 771, 771, 771, 772, 772, 772, 772, 772, 773, 773, 773, 773, 773, 774, 774, 774, 774, 774, 775, 775, 775, 775, 775, 776, 776, 776, 776, 776, 777, 777, 777, 777, 777, 778, 778, 778, 778, 778, 779, 779, 779, 779, 779, 780, 780, 780, 780, 780, 781, 781, 781, 781, 781, 782, 782, 782, 782, 782, 783, 783, 783, 783, 783, 784, 784, 784, 784, 784, 785, 785, 785, 785, 785, 786, 786, 786, 786, 786, 787, 787, 787, 787, 787, 788, 788, 788, 788, 788, 789, 789, 789, 789, 789, 790, 790, 790, 790, 790, 791, 791, 791, 791, 791, 792, 792, 792, 792, 792, 793, 793, 793, 793, 793, 794, 794, 794, 794, 794, 795, 795, 795, 795, 795, 796, 796, 796, 796, 796, 797, 797, 797, 797, 797, 798, 798, 798, 798, 798, 799, 799, 799, 799, 799, 800, 800, 800, 800, 800, 801, 801, 801, 801, 801, 802, 802, 802, 802, 802, 803, 803, 803, 803, 803, 804, 804, 804, 804, 804, 805, 805, 805, 805, 805, 806, 806, 806, 806, 806, 807, 807, 807, 807, 807, 808, 808, 808, 808, 808, 809, 809, 809, 809, 809, 810, 810, 810, 810, 810, 811, 811, 811, 811, 811, 812, 812, 812, 812, 812, 813, 813, 813, 813, 813, 814, 814, 814, 814, 814, 815, 815, 815, 815, 815, 816, 816, 816, 816, 816, 817, 817, 817, 817, 817, 818, 818, 818, 818, 818, 819, 819, 819, 819, 819, 820, 820, 820, 820, 820, 821, 821, 821, 821, 821, 822, 822, 822, 822, 822, 823, 823, 823, 823, 823, 824, 824, 824, 824, 824, 825, 825, 825, 825, 825, 826, 826, 826, 826, 826, 827, 827, 827, 827, 827, 828, 828, 828, 828, 828, 829, 829, 829, 829, 829, 830, 830, 830, 830, 830, 831, 831, 831, 831, 831, 832, 832, 832, 832, 832, 833, 833, 833, 833, 833, 834, 834, 834, 834, 834, 835, 835, 835, 835, 835, 836, 836, 836, 836, 836, 837, 837, 837, 837, 837, 838, 838, 838, 838, 838, 839, 839, 839, 839, 839, 840, 840, 840, 840, 840, 841, 841, 841, 841, 841, 842, 842, 842, 842, 842, 843, 843, 843, 843, 843, 844, 844, 844, 844, 844, 845, 845, 845, 845, 845, 846, 846, 846, 846, 846, 847, 847, 847, 847, 847, 848, 848, 848, 848, 848, 849, 849, 849, 849, 849, 850, 850, 850, 850, 850, 851, 851, 851, 851, 851, 852, 852, 852, 852, 852, 853, 853, 853, 853, 853, 854, 854, 854, 854, 854, 855, 855, 855, 855, 855, 856, 856, 856, 856, 856, 857, 857, 857, 857, 857, 858, 858, 858, 858, 858, 859, 859, 859, 859, 859, 860, 860, 860, 860, 860, 861, 861, 861, 861, 861, 862, 862, 862, 862, 862, 863, 863, 863, 863, 863, 864, 864, 864, 864, 864, 865, 865, 865, 865, 865, 866, 866, 866, 866, 866, 867, 867, 867, 867, 867, 868, 868, 868, 868, 868, 869, 869, 869, 869, 869, 870, 870, 870, 870, 870, 871, 871, 871, 871, 871, 872, 872, 872, 872, 872, 873, 873, 873, 873, 873, 874, 874, 874, 874, 874, 875, 875, 875, 875, 875, 876, 876, 876, 876, 876, 877, 877, 877, 877, 877, 878, 878, 878, 878, 878, 879, 879, 879, 879, 879, 880, 880, 880, 880, 880, 881, 881, 881, 881, 881, 882, 882, 882, 882, 882, 883, 883, 883, 883, 883, 884, 884, 884, 884, 884, 885, 885, 885, 885, 885, 886, 886, 886, 886, 886, 887, 887, 887, 887, 887, 888, 888, 888, 888, 888, 889, 889, 889, 889, 889, 890, 890, 890, 890, 890, 891, 891, 891, 891, 891, 892, 892, 892, 892, 892, 893, 893, 893, 893, 893, 894, 894, 894, 894, 894, 895, 895, 895, 895, 895, 896, 896, 896, 896, 896, 897, 897, 897, 897, 897, 898, 898, 898, 898, 898, 899, 899, 899, 899, 899, 900, 900, 900, 900, 900, 901, 901, 901, 901, 901, 902, 902, 902, 902, 902, 903, 903, 903, 903, 903, 904, 904, 904, 904, 904, 905, 905, 905, 905, 905, 906, 906, 906, 906, 906, 907, 907, 907, 907, 907, 908, 908, 908, 908, 908, 909, 909, 909, 909, 909, 910, 910, 910, 910, 910, 911, 911, 911, 911, 911, 912, 912, 912, 912, 912, 913, 913, 913, 913, 913, 914, 914, 914, 914, 914, 915, 915, 915, 915, 915, 916, 916, 916, 916, 916, 917, 917, 917, 917, 917, 918, 918, 918, 918, 918, 919, 919, 919, 919, 919, 920, 920, 920, 920, 920, 921, 921, 921, 921, 921, 922, 922, 922, 922, 922, 923, 923, 923, 923, 923, 924, 924, 924, 924, 924, 925, 925, 925, 925, 925, 926, 926, 926, 926, 926, 927, 927, 927, 927, 927, 928, 928, 928, 928, 928, 929, 929, 929, 929, 929, 930, 930, 930, 930, 930, 931, 931, 931, 931, 931, 932, 932, 932, 932, 932, 933, 933, 933, 933, 933, 934, 934, 934, 934, 934, 935, 935, 935, 935, 935, 936, 936, 936, 936, 936, 937, 937, 937, 937, 937, 938, 938, 938, 938, 938, 939, 939, 939, 939, 939, 940, 940, 940, 940, 940, 941, 941, 941, 941, 941, 942, 942, 942, 942, 942, 943, 943, 943, 943, 943, 944, 944, 944, 944, 944, 945, 945, 945, 945, 945, 946, 946, 946, 946, 946, 947, 947, 947, 947, 947, 948, 948, 948, 948, 948, 949, 949, 949, 949, 949, 950, 950, 950, 950, 950, 951, 951, 951, 951, 951, 952, 952, 952, 952, 952, 953, 953, 953, 953, 953, 954, 954, 954, 954, 954, 955, 955, 955, 955, 955, 956, 956, 956, 956, 956, 957, 957, 957, 957, 957, 958, 958, 958, 958, 958, 959, 959, 959, 959, 959, 960, 960, 960, 960, 960, 961, 961, 961, 961, 961, 962, 962, 962, 962, 962, 963, 963, 963, 963, 963, 964, 964, 964, 964, 964, 965, 965, 965, 965, 965, 966, 966, 966, 966, 966, 967, 967, 967, 967, 967, 968, 968, 968, 968, 968, 969, 969, 969, 969, 969, 970, 970, 970, 970, 970, 971, 971, 971, 971, 971, 972, 972, 972, 972, 972, 973, 973, 973, 973, 973, 974, 974, 974, 974, 974, 975, 975, 975, 975, 975, 976, 976, 976, 976, 976, 977, 977, 977, 977, 977, 978, 978, 978, 978, 978, 979, 979, 979, 979, 979, 980, 980, 980, 980, 980, 981, 981, 981, 981, 981, 982, 982, 982, 982, 982, 983, 983, 983, 983, 983, 984, 984, 984, 984, 984, 985, 985, 985, 985, 985, 986, 986, 986, 986, 986, 987, 987, 987, 987, 987, 988, 988, 988, 988, 988, 989, 989, 989, 989, 989, 990, 990, 990, 990, 990, 991, 991, 991, 991, 991, 992, 992, 992, 992, 992, 993, 993, 993, 993, 993, 994, 994, 994, 994, 994, 995, 995, 995, 995, 995, 996, 996, 996, 996, 996, 997, 997, 997, 997, 997, 998, 998, 998, 998, 998, 999, 999, 999, 999, 999, 1000, 1000, 1000, 1000, 1000, 1001, 1001, 1001, 1001, 1001, 1002, 1002, 1002, 1002, 1002, 1003, 1003, 1003, 1003, 1003, 1004, 1004, 1004, 1004, 1004, 1005, 1005, 1005, 1005, 1005, 1006, 1006, 1006, 1006, 1006, 1007, 1007, 1007, 1007, 1007, 1008, 1008, 1008, 1008, 1008, 1009, 1009, 1009, 1009, 1009, 1010, 1010, 1010, 1010, 1010, 1011, 1011, 1011, 1011, 1011, 1012, 1012, 1012, 1012, 1012, 1013, 1013, 1013, 1013, 1013, 1014, 1014, 1014, 1014, 1014, 1015, 1015, 1015, 1015, 1015, 1016, 1016, 1016, 1016, 1016, 1017, 1017, 1017, 1017, 1017, 1018, 1018, 1018, 1018, 1018, 1019, 1019, 1019, 1019, 1019, 1020, 1020, 1020, 1020, 1020, 1021, 1021, 1021, 1021, 1021, 1022, 1022, 1022, 1022, 1022, 1023, 1023, 1023, 1023, 1023, 1024, 1024, 1024, 1024, 1024], [4.0, -1.0, -1.0, -1.0, -1.0, -1.0, 4.0, -1.0, -1.0, -1.0  …  -1.0, -1.0, -1.0, 4.0, -1.0, -1.0, -1.0, -1.0, -1.0, 4.0], 1024, 1024), [1.0 0.0 … 0.0 0.0; 0.0 1.0 … 0.0 0.0; … ; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 0.0 1.0], [0.0 0.03842943919353914 … 0.15224093497742697 0.03842943919353936; 0.03842943919353914 0.07685887838707828 … 0.1906703741709661 0.0768588783870785; … ; 0.15224093497742697 0.1906703741709661 … 0.30448186995485393 0.19067037417096633; 0.03842943919353914 0.07685887838707828 … 0.1906703741709661 0.0768588783870785])

One of the benefits of the Bayesian approach is that we can fit for the hyperparameters of our prior/regularizers unlike traditional RML appraoches. To construct this heirarchical prior we will first make a map that takes in our regularizer hyperparameters and returns the image prior given those hyperparameters.

fmap = let crcache=crcache
    x->GaussMarkovRandomField(x, 1.0, crcache)
end

#1 (generic function with 1 method)

Now we can finally form our image prior. For this we use a heirarchical prior where the correlation length is given by a inverse gamma prior to prevent overfitting. Gaussian Markov random fields are extremly flexible models. To prevent overfitting it is common to use priors that penalize complexity. Therefore, we want to use priors that enforce similarity to our mean image, and prefer smoothness.

cprior = HierarchicalPrior(fmap, InverseGamma(1.0, -log(0.01*rat)))

prior = NamedDist(c = cprior, σimg = truncated(Normal(0.0, 1.0); lower=0.01), fg=Uniform(0.0, 1.0))

lklhd = RadioLikelihood(sky, dlcamp, dcphase;
                        skymeta = metadata)
post = Posterior(lklhd, prior)

Posterior{RadioLikelihood{Comrade.ModelMetadata{typeof(Main.sky), NamedTuple{(:meanpr, :K, :grid, :cache), Tuple{Matrix{Float64}, VLBIImagePriors.CenterImage{Matrix{Float64}, Tuple{Int64, Int64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}, VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}}}}, Nothing, Tuple{Comrade.EHTObservation{Float64, Comrade.EHTLogClosureAmplitudeDatum{Float64}, StructArrays.StructVector{Comrade.EHTLogClosureAmplitudeDatum{Float64}, NamedTuple{(:measurement, :error, :U1, :V1, :U2, :V2, :U3, :V3, :U4, :V4, :T, :F, :quadrangle), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{NTuple{4, Symbol}}}}, Int64}, Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}, Int64}, Comrade.EHTObservation{Float64, Comrade.EHTClosurePhaseDatum{Float64}, StructArrays.StructVector{Comrade.EHTClosurePhaseDatum{Float64}, NamedTuple{(:measurement, :error, :U1, :V1, :U2, :V2, :U3, :V3, :T, :F, :triangle), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol, Symbol}}}}, Int64}, Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}, Int64}}, Tuple{Comrade.ConditionedLikelihood{Comrade.var"#34#35"{Float64, Base.Fix2{typeof(logclosure_amplitudes), Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}}, VLBILikelihoods.CholeskyFactor{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.CHOLMOD.Factor{Float64}}}, Vector{Float64}}, Comrade.ConditionedLikelihood{Comrade.var"#36#37"{Float64, Base.Fix2{typeof(closure_phases), Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}}, VLBILikelihoods.CholeskyFactor{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.CHOLMOD.Factor{Float64}}}, Vector{Float64}}}, Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}, NamedTuple{(:U, :V, :T, :F), NTuple{4, Vector{Float64}}}}, VLBIImagePriors.NamedDist{(:c, :σimg, :fg), Tuple{VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Nothing}, Distributions.Uniform{Float64}}}}(RadioLikelihood
	Number of data products: 2
, VLBIImagePriors.NamedDist{(:c, :σimg, :fg), Tuple{VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Nothing}, Distributions.Uniform{Float64}}}(
dists: (VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}(
priormap: #1
hyperprior: Distributions.InverseGamma{Float64}(
invd: Distributions.Gamma{Float64}(α=1.0, θ=0.3410052124381751)
θ: 2.932506494109095
)

)
, Truncated(Distributions.Normal{Float64}(μ=0.0, σ=1.0); lower=0.01), Distributions.Uniform{Float64}(a=0.0, b=1.0))
)
)

Reconstructing the Image

To sample from this posterior, it is convenient to first move from our constrained parameter space to an unconstrained one (i.e., the support of the transformed posterior is (-∞, ∞)). This is done using the asflat function.

tpost = asflat(post)

Comrade.TransformedPosterior{Posterior{RadioLikelihood{Comrade.ModelMetadata{typeof(Main.sky), NamedTuple{(:meanpr, :K, :grid, :cache), Tuple{Matrix{Float64}, VLBIImagePriors.CenterImage{Matrix{Float64}, Tuple{Int64, Int64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}, VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}}}}, Nothing, Tuple{Comrade.EHTObservation{Float64, Comrade.EHTLogClosureAmplitudeDatum{Float64}, StructArrays.StructVector{Comrade.EHTLogClosureAmplitudeDatum{Float64}, NamedTuple{(:measurement, :error, :U1, :V1, :U2, :V2, :U3, :V3, :U4, :V4, :T, :F, :quadrangle), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{NTuple{4, Symbol}}}}, Int64}, Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}, Int64}, Comrade.EHTObservation{Float64, Comrade.EHTClosurePhaseDatum{Float64}, StructArrays.StructVector{Comrade.EHTClosurePhaseDatum{Float64}, NamedTuple{(:measurement, :error, :U1, :V1, :U2, :V2, :U3, :V3, :T, :F, :triangle), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol, Symbol}}}}, Int64}, Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}, Int64}}, Tuple{Comrade.ConditionedLikelihood{Comrade.var"#34#35"{Float64, Base.Fix2{typeof(logclosure_amplitudes), Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}}, VLBILikelihoods.CholeskyFactor{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.CHOLMOD.Factor{Float64}}}, Vector{Float64}}, Comrade.ConditionedLikelihood{Comrade.var"#36#37"{Float64, Base.Fix2{typeof(closure_phases), Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}}, VLBILikelihoods.CholeskyFactor{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.CHOLMOD.Factor{Float64}}}, Vector{Float64}}}, Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}, NamedTuple{(:U, :V, :T, :F), NTuple{4, Vector{Float64}}}}, VLBIImagePriors.NamedDist{(:c, :σimg, :fg), Tuple{VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Nothing}, Distributions.Uniform{Float64}}}}, TransformVariables.TransformTuple{NamedTuple{(:c, :σimg, :fg), Tuple{TransformVariables.TransformTuple{NamedTuple{(:params, :hyperparams), Tuple{TransformVariables.ArrayTransformation{TransformVariables.Identity, 2}, TransformVariables.ShiftedExp{true, Float64}}}}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ScaledShiftedLogistic{Float64}}}}}(Posterior{RadioLikelihood{Comrade.ModelMetadata{typeof(Main.sky), NamedTuple{(:meanpr, :K, :grid, :cache), Tuple{Matrix{Float64}, VLBIImagePriors.CenterImage{Matrix{Float64}, Tuple{Int64, Int64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}, VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}}}}, Nothing, Tuple{Comrade.EHTObservation{Float64, Comrade.EHTLogClosureAmplitudeDatum{Float64}, StructArrays.StructVector{Comrade.EHTLogClosureAmplitudeDatum{Float64}, NamedTuple{(:measurement, :error, :U1, :V1, :U2, :V2, :U3, :V3, :U4, :V4, :T, :F, :quadrangle), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{NTuple{4, Symbol}}}}, Int64}, Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}, Int64}, Comrade.EHTObservation{Float64, Comrade.EHTClosurePhaseDatum{Float64}, StructArrays.StructVector{Comrade.EHTClosurePhaseDatum{Float64}, NamedTuple{(:measurement, :error, :U1, :V1, :U2, :V2, :U3, :V3, :T, :F, :triangle), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol, Symbol}}}}, Int64}, Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}, Int64}}, Tuple{Comrade.ConditionedLikelihood{Comrade.var"#34#35"{Float64, Base.Fix2{typeof(logclosure_amplitudes), Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}}, VLBILikelihoods.CholeskyFactor{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.CHOLMOD.Factor{Float64}}}, Vector{Float64}}, Comrade.ConditionedLikelihood{Comrade.var"#36#37"{Float64, Base.Fix2{typeof(closure_phases), Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}}, VLBILikelihoods.CholeskyFactor{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.CHOLMOD.Factor{Float64}}}, Vector{Float64}}}, Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}, NamedTuple{(:U, :V, :T, :F), NTuple{4, Vector{Float64}}}}, VLBIImagePriors.NamedDist{(:c, :σimg, :fg), Tuple{VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Nothing}, Distributions.Uniform{Float64}}}}(RadioLikelihood
	Number of data products: 2
, VLBIImagePriors.NamedDist{(:c, :σimg, :fg), Tuple{VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Nothing}, Distributions.Uniform{Float64}}}(
dists: (VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}(
priormap: #1
hyperprior: Distributions.InverseGamma{Float64}(
invd: Distributions.Gamma{Float64}(α=1.0, θ=0.3410052124381751)
θ: 2.932506494109095
)

)
, Truncated(Distributions.Normal{Float64}(μ=0.0, σ=1.0); lower=0.01), Distributions.Uniform{Float64}(a=0.0, b=1.0))
)
), TransformVariables.TransformTuple{NamedTuple{(:c, :σimg, :fg), Tuple{TransformVariables.TransformTuple{NamedTuple{(:params, :hyperparams), Tuple{TransformVariables.ArrayTransformation{TransformVariables.Identity, 2}, TransformVariables.ShiftedExp{true, Float64}}}}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ScaledShiftedLogistic{Float64}}}}((c = TransformVariables.TransformTuple{NamedTuple{(:params, :hyperparams), Tuple{TransformVariables.ArrayTransformation{TransformVariables.Identity, 2}, TransformVariables.ShiftedExp{true, Float64}}}}((params = TransformVariables.ArrayTransformation{TransformVariables.Identity, 2}(asℝ, (32, 32)), hyperparams = as(Real, 0.0, ∞)), 1025), σimg = as(Real, 0.01, ∞), fg = as(Real, 0.0, 1.0)), 1027))

We can now also find the dimension of our posterior or the number of parameters we will sample.

ndim = dimension(tpost)

1027

Now we optimize using LBFGS

using ComradeOptimization
using OptimizationOptimJL
using Zygote
f = OptimizationFunction(tpost, Optimization.AutoZygote())
prob = Optimization.OptimizationProblem(f, prior_sample(rng, tpost), nothing)
sol = solve(prob, LBFGS(); maxiters=5_00);

Before we analyze our solution we first need to transform back to parameter space.

xopt = transform(tpost, sol)

(c = (params = [-0.8494620261149127 -0.6903432132227583 … -1.1460429104665082 -0.9962305690916786; -0.804566006281114 -0.634332548000261 … -1.283397291236575 -1.025321655049105; … ; -0.7968146139720274 -0.6670668747733154 … -0.9086999181453329 -0.8422587408717032; -0.8009085653610114 -0.6496714575767522 … -1.025180006630898 -0.9163586462477289], hyperparams = 19.018518327712204), σimg = 3.381626657616295, fg = 0.6883212825680473)

First we will evaluate our fit by plotting the residuals

using Plots
residual(skymodel(post, xopt), dlcamp, ylabel="Log Closure Amplitude Res.")

and now closure phases

residual(skymodel(post, xopt), dcphase, ylabel="|Closure Phase Res.|")

Now let's plot the MAP estimate.

import WGLMakie as CM
img = intensitymap(skymodel(post, xopt), μas2rad(150.0), μas2rad(150.0), 100, 100)
CM.image(img, axis=(xreversed=true, aspect=1, title="MAP Image"), colormap=:afmhot, figure=(;resolution=(400, 400),))

To sample from the posterior we will use HMC and more specifically the NUTS algorithm. For information about NUTS see Michael Betancourt's notes.

using ComradeAHMC
using Zygote
metric = DiagEuclideanMetric(ndim)
chain, stats = sample(post, AHMC(;metric, autodiff=Val(:Zygote)), 700; nadapts=500, init_params=xopt)

(NamedTuple{(:c, :σimg, :fg), Tuple{NamedTuple{(:params, :hyperparams), Tuple{Matrix{Float64}, Float64}}, Float64, Float64}}[(c = (params = [-0.9462351771077331 -0.6975258664214178 … -1.1730488151042615 -0.9449714870358811; -0.7816151663792513 -0.6407229074028492 … -1.3506955920070498 -0.9467244787760345; … ; -0.7489769331050818 -0.6441778166434293 … -0.8975879826019837 -0.9176535201985786; -0.7737310566097165 -0.6333205468635623 … -1.0087926598289623 -0.8730363780347268], hyperparams = 19.63619547827477), σimg = 3.494541487272257, fg = 0.6926728712080034), (c = (params = [-1.062602295274311 -0.18389402297487253 … -0.3386383289846316 -0.4155995520874266; -1.678745155800979 -0.35676851515820796 … -1.2988249531740057 -1.5039115585085114; … ; 0.7549584545763055 0.06353538431614432 … -0.5376905170741717 0.2823331683787087; -0.7422585871564134 -1.3905389295783637 … -0.8146345723832336 0.057053443176966354], hyperparams = 31.433891010036312), σimg = 2.1997348600378057, fg = 0.6309006376052073), (c = (params = [0.2812715016552185 -0.5539572695879331 … -1.2258487689604667 -0.8840335560587025; 0.62440144192357 -0.18635402122274133 … -0.4124736028475872 -0.08381871736477267; … ; -2.214924760467808 -1.678849110673396 … -0.5398643492488352 -1.9518154555890765; -0.23122911721742487 0.48181583301542635 … -1.0131695548664728 -0.7285131806649956], hyperparams = 38.04112327917421), σimg = 1.4149257130828912, fg = 0.5018565346983168), (c = (params = [-0.37006334607382757 0.07589001913729276 … -0.36951967036442224 0.06377385370615883; -0.11261735184379339 -0.20683181164891096 … -0.6918506812287294 0.4212916609314396; … ; -0.9740281328720336 0.2861797177809091 … -0.6457868599644515 -1.6931635949599457; -0.4037985422102937 -1.2552317623379436 … -0.9928537924629096 -1.8541859234816496], hyperparams = 12.615921829124828), σimg = 1.7768591698585163, fg = 0.5691240838167585), (c = (params = [-0.8405708525422453 -1.2740293752263836 … -0.9125766427169607 -0.742266787355985; -0.8629609548462537 -0.7118666738072563 … -0.16444408872810176 -0.4793999064911864; … ; -0.5834855577349618 -0.8460351944969285 … -0.8221980917053939 -0.27572318958181186; -1.1350996649248917 0.46996423572804774 … -1.0894287991143086 -0.0828841077700087], hyperparams = 106.59029468628931), σimg = 1.5546822227498114, fg = 0.3199936444422504), (c = (params = [-0.8435739616227437 -1.2860056089667373 … -0.9232481289760317 -0.7420965688212912; -0.8612886652929286 -0.7110236059811342 … -0.15735491135772664 -0.4809794847286285; … ; -0.5626689451200156 -0.85230966836672 … -0.8401004765716991 -0.27540922233789145; -1.156501036191873 0.4781241909889329 … -1.0794318058729646 -0.08125868032010827], hyperparams = 107.53394399079035), σimg = 1.5471324434266416, fg = 0.3192209947160163), (c = (params = [-0.7736903714655502 -0.8748262184852462 … -0.0709770905241182 -1.1228833487899028; -0.6261088011411556 -1.3915563512411357 … -0.8797026145362267 -0.609946605113679; … ; -1.011754712991703 0.14014947901134414 … -2.294940205429686 -1.8526673516607304; -0.026051770996812072 -1.6913963231035478 … -0.44377475811685596 -1.192011886084195], hyperparams = 57.059044873347894), σimg = 1.4026805079248506, fg = 0.4508014653102825), (c = (params = [-0.8801850710691077 -1.0331611080320786 … -0.5956009024204342 -0.6597392968190569; -0.7756684145529806 -1.1666980931429352 … -0.6331432599175869 -0.5769706972148662; … ; -1.2511218912540245 -0.09394595275457243 … -2.074272076620974 -1.9228192263793513; 0.18634957548526512 -1.212502794672401 … -1.08972399656909 -1.0125051746537141], hyperparams = 88.4613591717273), σimg = 1.813071740875319, fg = 0.5107236364206656), (c = (params = [-1.2365013215556486 -1.5995068431300956 … -1.0051498854898686 -1.067953697360107; -1.1681986670469586 -0.8109018209748344 … -0.8357215241123835 -1.062377955222338; … ; -0.3352996052956795 -0.23738056180094916 … -2.225144854757743 -1.4709545825903656; -1.7634835302865617 -0.7087851152151795 … -1.9099213006216447 -1.2609478153242661], hyperparams = 15.512515138059593), σimg = 1.8091654405959747, fg = 0.4852486974931797), (c = (params = [-1.2365013215556486 -1.5995068431300956 … -1.0051498854898686 -1.067953697360107; -1.1681986670469586 -0.8109018209748344 … -0.8357215241123835 -1.062377955222338; … ; -0.3352996052956795 -0.23738056180094916 … -2.225144854757743 -1.4709545825903656; -1.7634835302865617 -0.7087851152151795 … -1.9099213006216447 -1.2609478153242661], hyperparams = 15.512515138059593), σimg = 1.8091654405959747, fg = 0.4852486974931797)  …  (c = (params = [0.021335566742223533 0.17794721664477375 … 0.194311886638769 -0.3157697757338721; -0.836220523343804 -0.07559949551382433 … -0.8886612757597931 -0.06400137235845674; … ; 0.3773842578887056 0.6288240580057383 … 0.4995458396322825 -0.33854990401453205; -0.059802030203200116 1.4806150998386711 … 0.45971417646385404 -0.1863397294559373], hyperparams = 3.0486554869187454), σimg = 2.5070588544699324, fg = 0.3539302423839184), (c = (params = [-0.8620398231896101 -0.43234643975366027 … -0.21904293949392928 0.08058409481413344; -0.5215009512292117 -1.058999213952316 … 0.38020875811599625 -1.1535429545592095; … ; 0.01765388804188612 0.4219039517669606 … 0.014745657037428818 0.8077571347264895; 0.017382042498417374 -1.1747958975301345 … -0.05181499310763299 0.24951853738604515], hyperparams = 2.898371924934943), σimg = 2.380141752482626, fg = 0.3519783085779211), (c = (params = [-0.8243939527813051 -1.0411494009977411 … 0.09616187808570598 -1.1446848194218289; -1.6055802022476977 -0.6180997958672081 … 0.08653012984667749 -0.5787422077952221; … ; 0.20160460668998564 0.2656516140557881 … -0.059235621295065384 0.07472664402820574; -0.5713603237633886 0.3884208726742153 … 0.4370295088250859 0.4372002135814307], hyperparams = 4.715566519021752), σimg = 2.0632015639325663, fg = 0.5125859190125202), (c = (params = [-0.5688438445873175 0.4588629897433608 … 0.5934060353902498 0.12461662725701868; 0.003197651793439081 -0.486637426130007 … 0.08264607059618877 -0.7290632333923197; … ; 0.41577019919934954 0.535813112955021 … 0.6603561139908398 0.9816811140209425; 0.08184522206090301 -0.2927264462904561 … -0.25541518209210556 -0.60103082707821], hyperparams = 5.338486775398733), σimg = 1.935561996229206, fg = 0.4490710650915971), (c = (params = [-0.11078147898647932 -1.1675668022911603 … -0.7290454940157883 -0.6117491625898088; -0.460178389327492 -0.41089132182802446 … 0.639028443728355 0.8256857254263639; … ; -0.33236214246538937 -0.9934716388711774 … 0.678634561115883 0.1659715848636883; -1.00970640699836 -0.021975834864290354 … 0.22005083381887028 -0.3636238640605552], hyperparams = 5.304746012940809), σimg = 2.510383423018713, fg = 0.3842967194252401), (c = (params = [0.2933745941597561 0.7205806266222277 … -0.4003298965113413 0.5808461233967768; 0.6294313703532565 -0.8843809288358137 … 0.11840009711297422 0.2757050288267131; … ; 0.662387175922236 0.5566277112522631 … 0.6421580520643807 -0.025610965981084317; -0.04639780864979078 -0.10654780489723034 … -0.07544108237955054 0.3387769592914475], hyperparams = 5.84950336593517), σimg = 2.1458911774395952, fg = 0.44268173772118347), (c = (params = [-0.013933776767618724 -0.5447082400893417 … 0.4651493146354862 -1.205167107701235; 0.2877439699642336 1.4765459806249233 … 0.3682540152497222 0.0827961951669795; … ; -0.05352861851121698 0.08945538605977275 … -0.12461092398302426 -0.10616046615713967; 0.1901797693715514 0.14691991031767987 … 0.24255794298994637 -0.8746176404271078], hyperparams = 6.072999553809627), σimg = 2.2599863040879624, fg = 0.3080731586151542), (c = (params = [0.12730952532315462 -0.4569085907216185 … 0.3588876934782149 -1.089122752925769; 0.09624718564847075 1.3434802240508072 … 0.21798098448777853 0.1255799921264312; … ; -0.13831899248334012 0.11473604614419289 … -0.03994154287339072 -0.2083270684807581; 0.20636686955492275 -0.2928803992909045 … 0.26785633414707566 -0.9887785681106557], hyperparams = 5.151098670627499), σimg = 2.24262075548013, fg = 0.30810897898968814), (c = (params = [0.29428311437288174 -0.28744919743009145 … 0.13024875494215107 -1.204623606931678; 0.4005807435139919 1.4956995747409874 … -0.1211627463034794 0.34692875000419066; … ; 0.0031370950244578233 -0.03780338757106432 … 0.08403147247593652 -0.21928475120150054; 0.11316559280399835 -0.8504239802327918 … 0.15919638258148924 -0.7800312586774694], hyperparams = 3.57920601175801), σimg = 2.1882988208290883, fg = 0.3217049077767638), (c = (params = [-0.056101170027320166 0.658627444285125 … -0.6896382534519805 0.9079926532767154; -0.1627427521273243 -0.4324886488016908 … -0.1243276071929269 -0.4992945005338458; … ; 0.24072361241022122 0.2534432177265726 … -0.4934970789607874 -0.19589466626108187; -0.1731996507106475 1.274003674845524 … -0.8889983139435934 0.4332679112447181], hyperparams = 6.031582019850623), σimg = 1.712524099391445, fg = 0.09998331865290522)], NamedTuple{(:n_steps, :is_accept, :acceptance_rate, :log_density, :hamiltonian_energy, :hamiltonian_energy_error, :max_hamiltonian_energy_error, :tree_depth, :numerical_error, :step_size, :nom_step_size, :is_adapt), Tuple{Int64, Bool, Float64, Float64, Float64, Float64, Float64, Int64, Bool, Float64, Float64, Bool}}[(n_steps = 1023, is_accept = 1, acceptance_rate = 0.9992902834360256, log_density = -471.28034446776434, hamiltonian_energy = 963.0210610573565, hamiltonian_energy_error = 0.000802326070242998, max_hamiltonian_energy_error = 0.0015327301232446189, tree_depth = 10, numerical_error = 0, step_size = 0.0001, nom_step_size = 0.0001, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9251029951353055, log_density = -658.3531507868277, hamiltonian_energy = 996.5378095988806, hamiltonian_energy_error = 0.014666115759155218, max_hamiltonian_energy_error = -0.5345113532648611, tree_depth = 10, numerical_error = 0, step_size = 0.0015734254544819278, nom_step_size = 0.0015734254544819278, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.6137859757560344, log_density = -700.2500982979964, hamiltonian_energy = 1193.807827297034, hamiltonian_energy_error = 0.5879367981322048, max_hamiltonian_energy_error = 1.2575830648770534, tree_depth = 10, numerical_error = 0, step_size = 0.0027191098087846807, nom_step_size = 0.0027191098087846807, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9843949456582975, log_density = -671.9100847577166, hamiltonian_energy = 1222.7838611874058, hamiltonian_energy_error = -0.010088161882777058, max_hamiltonian_energy_error = -0.08894799292966127, tree_depth = 10, numerical_error = 0, step_size = 0.0021552557635181584, nom_step_size = 0.0021552557635181584, is_adapt = 1), (n_steps = 511, is_accept = 1, acceptance_rate = 0.9917889526755399, log_density = -651.720211612358, hamiltonian_energy = 1194.8634173848823, hamiltonian_energy_error = -0.32801178340696424, max_hamiltonian_energy_error = -0.43245519085394335, tree_depth = 9, numerical_error = 0, step_size = 0.00445073451815245, nom_step_size = 0.00445073451815245, is_adapt = 1), (n_steps = 6, is_accept = 1, acceptance_rate = 0.19379703186343591, log_density = -653.8505262498163, hamiltonian_energy = 1204.430214480151, hamiltonian_energy_error = -0.1327353351259717, max_hamiltonian_energy_error = 4713.731139979011, tree_depth = 2, numerical_error = 1, step_size = 0.009765840419875316, nom_step_size = 0.009765840419875316, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9952006433481879, log_density = -654.0012374069169, hamiltonian_energy = 1171.3378443313686, hamiltonian_energy_error = -0.004578211530997578, max_hamiltonian_energy_error = 0.01737291207950875, tree_depth = 10, numerical_error = 0, step_size = 0.0018914766041249179, nom_step_size = 0.0018914766041249179, is_adapt = 1), (n_steps = 511, is_accept = 1, acceptance_rate = 0.9374939431593299, log_density = -675.4424745670314, hamiltonian_energy = 1161.5024955125584, hamiltonian_energy_error = 0.07993217897160321, max_hamiltonian_energy_error = -0.35104153816132566, tree_depth = 9, numerical_error = 0, step_size = 0.004100680545900643, nom_step_size = 0.004100680545900643, is_adapt = 1), (n_steps = 255, is_accept = 1, acceptance_rate = 0.951548092630774, log_density = -612.865526781077, hamiltonian_energy = 1160.2813433825481, hamiltonian_energy_error = -0.31122515203151124, max_hamiltonian_energy_error = -0.9994830251566782, tree_depth = 8, numerical_error = 0, step_size = 0.007493322132379629, nom_step_size = 0.007493322132379629, is_adapt = 1), (n_steps = 3, is_accept = 1, acceptance_rate = 0.0004060093434629848, log_density = -612.865526781077, hamiltonian_energy = 1093.6591898554839, hamiltonian_energy_error = 0.0, max_hamiltonian_energy_error = 2134.8893560743313, tree_depth = 1, numerical_error = 1, step_size = 0.014299428471824607, nom_step_size = 0.014299428471824607, is_adapt = 1)  …  (n_steps = 255, is_accept = 1, acceptance_rate = 0.8128369575944936, log_density = -595.2358403071913, hamiltonian_energy = 1104.805252887043, hamiltonian_energy_error = -0.13887097842189178, max_hamiltonian_energy_error = 0.9320794108627979, tree_depth = 8, numerical_error = 0, step_size = 0.010772679071944955, nom_step_size = 0.010772679071944955, is_adapt = 0), (n_steps = 255, is_accept = 1, acceptance_rate = 0.8579360111943513, log_density = -600.8453029539152, hamiltonian_energy = 1103.374058984637, hamiltonian_energy_error = -0.5455914256529013, max_hamiltonian_energy_error = 1.6575360337049005, tree_depth = 8, numerical_error = 0, step_size = 0.010772679071944955, nom_step_size = 0.010772679071944955, is_adapt = 0), (n_steps = 255, is_accept = 1, acceptance_rate = 0.7561650052319309, log_density = -603.8413879111934, hamiltonian_energy = 1127.3427551292473, hamiltonian_energy_error = 0.22347481766951205, max_hamiltonian_energy_error = 1.2014543280688486, tree_depth = 8, numerical_error = 0, step_size = 0.010772679071944955, nom_step_size = 0.010772679071944955, is_adapt = 0), (n_steps = 255, is_accept = 1, acceptance_rate = 0.9028731825693388, log_density = -644.800145080576, hamiltonian_energy = 1111.2583141445466, hamiltonian_energy_error = -0.021773877862869995, max_hamiltonian_energy_error = -0.6195969656237139, tree_depth = 8, numerical_error = 0, step_size = 0.010772679071944955, nom_step_size = 0.010772679071944955, is_adapt = 0), (n_steps = 255, is_accept = 1, acceptance_rate = 0.890977507397043, log_density = -611.6194112069845, hamiltonian_energy = 1152.1960974881097, hamiltonian_energy_error = -0.15656378508424496, max_hamiltonian_energy_error = 2.9567809531195053, tree_depth = 8, numerical_error = 0, step_size = 0.010772679071944955, nom_step_size = 0.010772679071944955, is_adapt = 0), (n_steps = 255, is_accept = 1, acceptance_rate = 0.9093755378732006, log_density = -632.2930145917305, hamiltonian_energy = 1126.9719432387024, hamiltonian_energy_error = -0.21379628620729818, max_hamiltonian_energy_error = -0.6191212249268574, tree_depth = 8, numerical_error = 0, step_size = 0.010772679071944955, nom_step_size = 0.010772679071944955, is_adapt = 0), (n_steps = 255, is_accept = 1, acceptance_rate = 0.8593128880617118, log_density = -634.785078220036, hamiltonian_energy = 1155.0171384370356, hamiltonian_energy_error = 0.4666349606463882, max_hamiltonian_energy_error = 0.7416166697075823, tree_depth = 8, numerical_error = 0, step_size = 0.010772679071944955, nom_step_size = 0.010772679071944955, is_adapt = 0), (n_steps = 255, is_accept = 1, acceptance_rate = 0.6007128239617269, log_density = -635.1073248038476, hamiltonian_energy = 1154.1266570107457, hamiltonian_energy_error = -0.7421072576626102, max_hamiltonian_energy_error = 1.6378420275475492, tree_depth = 8, numerical_error = 0, step_size = 0.010772679071944955, nom_step_size = 0.010772679071944955, is_adapt = 0), (n_steps = 255, is_accept = 1, acceptance_rate = 0.8433004700697071, log_density = -610.8903435164816, hamiltonian_energy = 1124.8393153700308, hamiltonian_energy_error = 0.28282435304731735, max_hamiltonian_energy_error = 0.8837121971480428, tree_depth = 8, numerical_error = 0, step_size = 0.010772679071944955, nom_step_size = 0.010772679071944955, is_adapt = 0), (n_steps = 255, is_accept = 1, acceptance_rate = 0.7618050229729341, log_density = -603.2515767372013, hamiltonian_energy = 1102.7919372841984, hamiltonian_energy_error = 0.16317661893049262, max_hamiltonian_energy_error = 1.3317890751427512, tree_depth = 8, numerical_error = 0, step_size = 0.010772679071944955, nom_step_size = 0.010772679071944955, is_adapt = 0)])

Now that we have our posterior, we can assess which parts of the image are strongly inferred by the data. This is rather unique to Comrade where more traditional imaging algorithms like CLEAN and RML are inherently unable to assess uncertainty in their reconstructions.

To explore our posterior let's first create images from a bunch of draws from the posterior

msamples = skymodel.(Ref(post), chain[501:2:end]);

The mean image is then given by

using StatsBase
imgs = intensitymap.(msamples, μas2rad(150.0), μas2rad(150.0), 128, 128)
mimg = mean(imgs)
simg = std(imgs)
fig = CM.Figure(;resolution=(400, 400));
CM.image(fig[1,1], mimg,
                   axis=(xreversed=true, aspect=1, title="Mean Image"),
                   colormap=:afmhot)
CM.image(fig[1,2], simg./(max.(mimg, 1e-5)),
                   axis=(xreversed=true, aspect=1, title="1/SNR",), colorrange=(0.0, 2.0),
                   colormap=:afmhot)
CM.image(fig[2,1], imgs[1],
                   axis=(xreversed=true, aspect=1,title="Draw 1"),
                   colormap=:afmhot)
CM.image(fig[2,2], imgs[end],
                   axis=(xreversed=true, aspect=1,title="Draw 2"),
                   colormap=:afmhot)
fig

Now let's see whether our residuals look better.

p = plot();
for s in sample(chain[501:end], 10)
    residual!(p, vlbimodel(post, s), dlcamp)
end
ylabel!("Log-Closure Amplitude Res.");
p

p = plot();
for s in sample(chain[501:end], 10)
    residual!(p, vlbimodel(post, s), dcphase)
end
ylabel!("|Closure Phase Res.|");
p

And viola, you have a quick and preliminary image of M87 fitting only closure products. For a publication-level version we would recommend

1. Running the chain longer and multiple times to properly assess things like ESS and R̂ (see Geometric Modeling of EHT Data)
2. Fitting gains. Typically gain amplitudes are good to 10-20% for the EHT not the infinite uncertainty closures implicitly assume
3. Making sure the posterior is unimodal (hint for this example it isn't!). The EHT image posteriors can be pretty complicated, so typically you want to use a sampler that can deal with multi-modal posteriors. Check out the package Pigeons.jl for an in-development package that should easily enable this type of sampling.

This page was generated using Literate.jl.